home *** CD-ROM | disk | FTP | other *** search
/ AmigActive 21 / AACD 21.iso / AACD / Utilities / Ghostscript / src / gsmisc.c < prev    next >
Encoding:
C/C++ Source or Header  |  2001-01-01  |  30.1 KB  |  1,131 lines
  1. /* Copyright (C) 1989, 2000 Aladdin Enterprises.  All rights reserved.
  2.   
  3.   This file is part of AFPL Ghostscript.
  4.   
  5.   AFPL Ghostscript is distributed with NO WARRANTY OF ANY KIND.  No author or
  6.   distributor accepts any responsibility for the consequences of using it, or
  7.   for whether it serves any particular purpose or works at all, unless he or
  8.   she says so in writing.  Refer to the Aladdin Free Public License (the
  9.   "License") for full details.
  10.   
  11.   Every copy of AFPL Ghostscript must include a copy of the License, normally
  12.   in a plain ASCII text file named PUBLIC.  The License grants you the right
  13.   to copy, modify and redistribute AFPL Ghostscript, but only under certain
  14.   conditions described in the License.  Among other things, the License
  15.   requires that the copyright notice and this notice be preserved on all
  16.   copies.
  17. */
  18.  
  19. /*$Id: gsmisc.c,v 1.8 2000/09/19 19:00:30 lpd Exp $ */
  20. /* Miscellaneous utilities for Ghostscript library */
  21.  
  22. /*
  23.  * In order to capture the original definition of sqrt, which might be
  24.  * either a procedure or a macro and might not have an ANSI-compliant
  25.  * prototype (!), we need to do the following:
  26.  */
  27. #include "std.h"
  28. #if defined(VMS) && defined(__GNUC__)
  29. /*  DEC VAX/VMS C comes with a math.h file, but GNU VAX/VMS C does not. */
  30. #  include "vmsmath.h"
  31. #else
  32. #  include <math.h>
  33. #endif
  34. inline private double
  35. orig_sqrt(double x)
  36. {
  37.     return sqrt(x);
  38. }
  39.  
  40. /* Here is the real #include section. */
  41. #include "ctype_.h"
  42. #include "malloc_.h"
  43. #include "math_.h"
  44. #include "memory_.h"
  45. #include "string_.h"
  46. #include "gx.h"
  47. #include "gpcheck.h"        /* for gs_return_check_interrupt */
  48. #include "gserror.h"        /* for prototype */
  49. #include "gserrors.h"
  50. #include "gconfigv.h"        /* for USE_ASM */
  51. #include "gxfarith.h"
  52. #include "gxfixed.h"
  53.  
  54. /* Define private replacements for stdin, stdout, and stderr. */
  55. FILE *gs_stdio[3];
  56.  
  57. /* ------ Debugging ------ */
  58.  
  59. /* Ghostscript writes debugging output to gs_debug_out. */
  60. /* We define gs_debug and gs_debug_out even if DEBUG isn't defined, */
  61. /* so that we can compile individual modules with DEBUG set. */
  62. char gs_debug[128];
  63. FILE *gs_debug_out;
  64.  
  65. /* Test whether a given debugging option is selected. */
  66. /* Upper-case letters automatically include their lower-case counterpart. */
  67. bool
  68. gs_debug_c(int c)
  69. {
  70.     return
  71.     (c >= 'a' && c <= 'z' ? gs_debug[c] | gs_debug[c ^ 32] : gs_debug[c]);
  72. }
  73.  
  74. /* Define the formats for debugging printout. */
  75. const char *const dprintf_file_and_line_format = "%10s(%4d): ";
  76. const char *const dprintf_file_only_format = "%10s(unkn): ";
  77.  
  78. /*
  79.  * Define the trace printout procedures.  We always include these, in case
  80.  * other modules were compiled with DEBUG set.  Note that they must use
  81.  * fprintf, not fput[cs], because of the way that stdout is implemented on
  82.  * Windows platforms.
  83.  */
  84. void
  85. dflush(void)
  86. {
  87.     fflush(dstderr);
  88. }
  89. private const char *
  90. dprintf_file_tail(const char *file)
  91. {
  92.     const char *tail = file + strlen(file);
  93.  
  94.     while (tail > file &&
  95.        (isalnum(tail[-1]) || tail[-1] == '.' || tail[-1] == '_')
  96.     )
  97.     --tail;
  98.     return tail;
  99. }
  100. #if __LINE__            /* compiler provides it */
  101. void
  102. dprintf_file_and_line(FILE * f, const char *file, int line)
  103. {
  104.     if (gs_debug['/'])
  105.     fprintf(f, dprintf_file_and_line_format,
  106.         dprintf_file_tail(file), line);
  107. }
  108. #else
  109. void
  110. dprintf_file_only(FILE * f, const char *file)
  111. {
  112.     if (gs_debug['/'])
  113.     fprintf(f, dprintf_file_only_format, dprintf_file_tail(file));
  114. }
  115. #endif
  116. void
  117. printf_program_ident(FILE * f, const char *program_name,
  118.              long revision_number)
  119. {
  120.     if (program_name)
  121.     fprintf(f, (revision_number ? "%s " : "%s"), program_name);
  122.     if (revision_number) {
  123.     int fpart = revision_number % 100;
  124.  
  125.     fprintf(f, (fpart == 0 ? "%d.%d" : "%d.%02d"),
  126.         (int)(revision_number / 100), fpart);
  127.     }
  128. }
  129. void
  130. eprintf_program_ident(FILE * f, const char *program_name,
  131.               long revision_number)
  132. {
  133.     if (program_name) {
  134.     printf_program_ident(f, program_name, revision_number);
  135.     fprintf(f, ": ");
  136.     }
  137. }
  138. #if __LINE__            /* compiler provides it */
  139. void
  140. lprintf_file_and_line(FILE * f, const char *file, int line)
  141. {
  142.     fprintf(f, "%s(%d): ", file, line);
  143. }
  144. #else
  145. void
  146. lprintf_file_only(FILE * f, const char *file)
  147. {
  148.     fprintf(f, "%s(?): ", file);
  149. }
  150. #endif
  151.  
  152. /* Log an error return.  We always include this, in case other */
  153. /* modules were compiled with DEBUG set. */
  154. #undef gs_log_error        /* in case DEBUG isn't set */
  155. int
  156. gs_log_error(int err, const char *file, int line)
  157. {
  158.     if (gs_log_errors) {
  159.     if (file == NULL)
  160.         dprintf1("Returning error %d.\n", err);
  161.     else
  162.         dprintf3("%s(%d): Returning error %d.\n",
  163.              (const char *)file, line, err);
  164.     }
  165.     return err;
  166. }
  167.  
  168. /* Check for interrupts before a return. */
  169. int
  170. gs_return_check_interrupt(int code)
  171. {
  172.     if (code < 0)
  173.     return code;
  174.     {
  175.     int icode = gp_check_interrupts();
  176.  
  177.     return (icode == 0 ? code :
  178.         gs_note_error((icode > 0 ? gs_error_interrupt : icode)));
  179.     }
  180. }
  181.  
  182. /* ------ Substitutes for missing C library functions ------ */
  183.  
  184. #ifdef MEMORY__NEED_MEMMOVE    /* see memory_.h */
  185. /* Copy bytes like memcpy, guaranteed to handle overlap correctly. */
  186. /* ANSI C defines the returned value as being the src argument, */
  187. /* but with the const restriction removed! */
  188. void *
  189. gs_memmove(void *dest, const void *src, size_t len)
  190. {
  191.     if (!len)
  192.     return (void *)src;
  193. #define bdest ((byte *)dest)
  194. #define bsrc ((const byte *)src)
  195.     /* We use len-1 for comparisons because adding len */
  196.     /* might produce an offset overflow on segmented systems. */
  197.     if (PTR_LE(bdest, bsrc)) {
  198.     register byte *end = bdest + (len - 1);
  199.  
  200.     if (PTR_LE(bsrc, end)) {
  201.         /* Source overlaps destination from above. */
  202.         register const byte *from = bsrc;
  203.         register byte *to = bdest;
  204.  
  205.         for (;;) {
  206.         *to = *from;
  207.         if (to >= end)    /* faster than = */
  208.             return (void *)src;
  209.         to++;
  210.         from++;
  211.         }
  212.     }
  213.     } else {
  214.     register const byte *from = bsrc + (len - 1);
  215.  
  216.     if (PTR_LE(bdest, from)) {
  217.         /* Source overlaps destination from below. */
  218.         register const byte *end = bsrc;
  219.         register byte *to = bdest + (len - 1);
  220.  
  221.         for (;;) {
  222.         *to = *from;
  223.         if (from <= end)    /* faster than = */
  224.             return (void *)src;
  225.         to--;
  226.         from--;
  227.         }
  228.     }
  229.     }
  230. #undef bdest
  231. #undef bsrc
  232.     /* No overlap, it's safe to use memcpy. */
  233.     memcpy(dest, src, len);
  234.     return (void *)src;
  235. }
  236. #endif
  237.  
  238. #ifdef MEMORY__NEED_MEMCPY    /* see memory_.h */
  239. void *
  240. gs_memcpy(void *dest, const void *src, size_t len)
  241. {
  242.     if (len > 0) {
  243. #define bdest ((byte *)dest)
  244. #define bsrc ((const byte *)src)
  245.     /* We can optimize this much better later on. */
  246.     register byte *end = bdest + (len - 1);
  247.     register const byte *from = bsrc;
  248.     register byte *to = bdest;
  249.  
  250.     for (;;) {
  251.         *to = *from;
  252.         if (to >= end)    /* faster than = */
  253.         break;
  254.         to++;
  255.         from++;
  256.     }
  257.     }
  258. #undef bdest
  259. #undef bsrc
  260.     return (void *)src;
  261. }
  262. #endif
  263.  
  264. #ifdef MEMORY__NEED_MEMCHR    /* see memory_.h */
  265. /* ch should obviously be char rather than int, */
  266. /* but the ANSI standard declaration uses int. */
  267. void *
  268. gs_memchr(const void *ptr, int ch, size_t len)
  269. {
  270.     if (len > 0) {
  271.     register const char *p = ptr;
  272.     register uint count = len;
  273.  
  274.     do {
  275.         if (*p == (char)ch)
  276.         return (void *)p;
  277.         p++;
  278.     } while (--count);
  279.     }
  280.     return 0;
  281. }
  282. #endif
  283.  
  284. #ifdef MEMORY__NEED_MEMSET    /* see memory_.h */
  285. /* ch should obviously be char rather than int, */
  286. /* but the ANSI standard declaration uses int. */
  287. void *
  288. gs_memset(void *dest, register int ch, size_t len)
  289. {
  290.     /*
  291.      * This procedure is used a lot to fill large regions of images,
  292.      * so we take some trouble to optimize it.
  293.      */
  294.     register char *p = dest;
  295.     register size_t count = len;
  296.  
  297.     ch &= 255;
  298.     if (len >= sizeof(long) * 3) {
  299.     long wd = (ch << 24) | (ch << 16) | (ch << 8) | ch;
  300.  
  301.     while (ALIGNMENT_MOD(p, sizeof(long)))
  302.         *p++ = (char)ch, --count;
  303.     for (; count >= sizeof(long) * 4;
  304.          p += sizeof(long) * 4, count -= sizeof(long) * 4
  305.          )
  306.         ((long *)p)[3] = ((long *)p)[2] = ((long *)p)[1] =
  307.         ((long *)p)[0] = wd;
  308.     switch (count >> ARCH_LOG2_SIZEOF_LONG) {
  309.     case 3:
  310.         *((long *)p) = wd; p += sizeof(long);
  311.     case 2:
  312.         *((long *)p) = wd; p += sizeof(long);
  313.     case 1:
  314.         *((long *)p) = wd; p += sizeof(long);
  315.         count &= sizeof(long) - 1;
  316.     case 0:
  317.     default:        /* can't happen */
  318.         DO_NOTHING;
  319.     }
  320.     }
  321.     /* Do any leftover bytes. */
  322.     for (; count > 0; --count)
  323.     *p++ = (char)ch;
  324.     return dest;
  325. }
  326. #endif
  327.  
  328. #ifdef malloc__need_realloc    /* see malloc_.h */
  329. /* Some systems have non-working implementations of realloc. */
  330. void *
  331. gs_realloc(void *old_ptr, size_t old_size, size_t new_size)
  332. {
  333.     void *new_ptr;
  334.  
  335.     if (new_size) {
  336.     new_ptr = malloc(new_size);
  337.     if (new_ptr == NULL)
  338.         return NULL;
  339.     } else
  340.     new_ptr = NULL;
  341.     /* We have to pass in the old size, since we have no way to */
  342.     /* determine it otherwise. */
  343.     if (old_ptr != NULL) {
  344.     if (new_ptr != NULL)
  345.         memcpy(new_ptr, old_ptr, min(old_size, new_size));
  346.     free(old_ptr);
  347.     }
  348.     return new_ptr;
  349. }
  350. #endif
  351.  
  352. /* ------ Debugging support ------ */
  353.  
  354. /* Dump a region of memory. */
  355. void
  356. debug_dump_bytes(const byte * from, const byte * to, const char *msg)
  357. {
  358.     const byte *p = from;
  359.  
  360.     if (from < to && msg)
  361.     dprintf1("%s:\n", msg);
  362.     while (p != to) {
  363.     const byte *q = min(p + 16, to);
  364.  
  365.     dprintf1("0x%lx:", (ulong) p);
  366.     while (p != q)
  367.         dprintf1(" %02x", *p++);
  368.     dputc('\n');
  369.     }
  370. }
  371.  
  372. /* Dump a bitmap. */
  373. void
  374. debug_dump_bitmap(const byte * bits, uint raster, uint height, const char *msg)
  375. {
  376.     uint y;
  377.     const byte *data = bits;
  378.  
  379.     for (y = 0; y < height; ++y, data += raster)
  380.     debug_dump_bytes(data, data + raster, (y == 0 ? msg : NULL));
  381. }
  382.  
  383. /* Print a string. */
  384. void
  385. debug_print_string(const byte * chrs, uint len)
  386. {
  387.     uint i;
  388.  
  389.     for (i = 0; i < len; i++)
  390.     dputc(chrs[i]);
  391.     dflush();
  392. }
  393.  
  394. /*
  395.  * The following code prints a hex stack backtrace on Linux/Intel systems.
  396.  * It is here to be patched into places where we need to print such a trace
  397.  * because of gdb's inability to put breakpoints in dynamically created
  398.  * threads.
  399.  *
  400.  * first_arg is the first argument of the procedure into which this code
  401.  * is patched.
  402.  */
  403. #define BACKTRACE(first_arg)\
  404.   BEGIN\
  405.     ulong *fp_ = (ulong *)&first_arg - 2;\
  406.     for (; fp_ && (fp_[1] & 0xff000000) == 0x08000000; fp_ = (ulong *)*fp_)\
  407.     dprintf2("  fp=0x%lx ip=0x%lx\n", (ulong)fp_, fp_[1]);\
  408.   END
  409.  
  410. /* ------ Arithmetic ------ */
  411.  
  412. /* Compute M modulo N.  Requires N > 0; guarantees 0 <= imod(M,N) < N, */
  413. /* regardless of the whims of the % operator for negative operands. */
  414. int
  415. imod(int m, int n)
  416. {
  417.     if (n <= 0)
  418.     return 0;        /* sanity check */
  419.     if (m >= 0)
  420.     return m % n;
  421.     {
  422.     int r = -m % n;
  423.  
  424.     return (r == 0 ? 0 : n - r);
  425.     }
  426. }
  427.  
  428. /* Compute the GCD of two integers. */
  429. int
  430. igcd(int x, int y)
  431. {
  432.     int c = x, d = y;
  433.  
  434.     if (c < 0)
  435.     c = -c;
  436.     if (d < 0)
  437.     d = -d;
  438.     while (c != 0 && d != 0)
  439.     if (c > d)
  440.         c %= d;
  441.     else
  442.         d %= c;
  443.     return d + c;        /* at most one is non-zero */
  444. }
  445.  
  446. /* Compute X such that A*X = B mod M.  See gxarith.h for details. */
  447. int
  448. idivmod(int a, int b, int m)
  449. {
  450.     /*
  451.      * Use the approach indicated in Knuth vol. 2, section 4.5.2, Algorithm
  452.      * X (p. 302) and exercise 15 (p. 315, solution p. 523).
  453.      */
  454.     int u1 = 0, u3 = m;
  455.     int v1 = 1, v3 = a;
  456.     /*
  457.      * The following loop will terminate with a * u1 = gcd(a, m) mod m.
  458.      * Then x = u1 * b / gcd(a, m) mod m.  Since we require that
  459.      * gcd(a, m) | gcd(a, b), it follows that gcd(a, m) | b, so the
  460.      * division is exact.
  461.      */
  462.     while (v3) {
  463.     int q = u3 / v3, t;
  464.  
  465.     t = u1 - v1 * q, u1 = v1, v1 = t;
  466.     t = u3 - v3 * q, u3 = v3, v3 = t;
  467.     }
  468.     return imod(u1 * b / igcd(a, m), m);
  469. }
  470.  
  471. /* Compute floor(log2(N)).  Requires N > 0. */
  472. int
  473. ilog2(int n)
  474. {
  475.     int m = n, l = 0;
  476.  
  477.     while (m >= 16)
  478.     m >>= 4, l += 4;
  479.     return
  480.     (m <= 1 ? 0 :
  481.      "\000\000\001\001\002\002\002\002\003\003\003\003\003\003\003\003"[m] + l);
  482. }
  483.  
  484. #if defined(NEED_SET_FMUL2FIXED) && !USE_ASM
  485.  
  486. /*
  487.  * Floating multiply with fixed result, for avoiding floating point in
  488.  * common coordinate transformations.  Assumes IEEE representation,
  489.  * 16-bit short, 32-bit long.  Optimized for the case where the first
  490.  * operand has no more than 16 mantissa bits, e.g., where it is a user space
  491.  * coordinate (which are often integers).
  492.  *
  493.  * The assembly language version of this code is actually faster than
  494.  * the FPU, if the code is compiled with FPU_TYPE=0 (which requires taking
  495.  * a trap on every FPU operation).  If there is no FPU, the assembly
  496.  * language version of this code is over 10 times as fast as the emulated FPU.
  497.  */
  498. int
  499. set_fmul2fixed_(fixed * pr, long /*float */ a, long /*float */ b)
  500. {
  501.     ulong ma = (ushort)(a >> 8) | 0x8000;
  502.     ulong mb = (ushort)(b >> 8) | 0x8000;
  503.     int e = 260 + _fixed_shift - ( ((byte)(a >> 23)) + ((byte)(b >> 23)) );
  504.     ulong p1 = ma * (b & 0xff);
  505.     ulong p = ma * mb;
  506.  
  507. #define p_bits (size_of(p) * 8)
  508.  
  509.     if ((byte) a) {        /* >16 mantissa bits */
  510.     ulong p2 = (a & 0xff) * mb;
  511.  
  512.     p += ((((uint) (byte) a * (uint) (byte) b) >> 8) + p1 + p2) >> 8;
  513.     } else
  514.     p += p1 >> 8;
  515.     if ((uint) e < p_bits)    /* e = -1 is possible */
  516.     p >>= e;
  517.     else if (e >= p_bits) {    /* also detects a=0 or b=0 */
  518.     *pr = fixed_0;
  519.     return 0;
  520.     } else if (e >= -(p_bits - 1) || p >= 1L << (p_bits - 1 + e))
  521.     return_error(gs_error_limitcheck);
  522.     else
  523.     p <<= -e;
  524.     *pr = ((a ^ b) < 0 ? -p : p);
  525.     return 0;
  526. }
  527. int
  528. set_dfmul2fixed_(fixed * pr, ulong /*double lo */ xalo, long /*float */ b, long /*double hi */ xahi)
  529. {
  530.     return set_fmul2fixed_(pr,
  531.                (xahi & (3L << 30)) +
  532.                ((xahi << 3) & 0x3ffffff8) +
  533.                (xalo >> 29),
  534.                b);
  535. }
  536.  
  537. #endif /* NEED_SET_FMUL2FIXED */
  538.  
  539. #if USE_FPU_FIXED
  540.  
  541. /*
  542.  * Convert from floating point to fixed point with scaling.
  543.  * These are efficient algorithms for FPU-less machines.
  544.  */
  545. #define mbits_float 23
  546. #define mbits_double 20
  547. int
  548. set_float2fixed_(fixed * pr, long /*float */ vf, int frac_bits)
  549. {
  550.     fixed mantissa;
  551.     int shift;
  552.  
  553.     if (!(vf & 0x7f800000)) {
  554.     *pr = fixed_0;
  555.     return 0;
  556.     }
  557.     mantissa = (fixed) ((vf & 0x7fffff) | 0x800000);
  558.     shift = ((vf >> 23) & 255) - (127 + 23) + frac_bits;
  559.     if (shift >= 0) {
  560.     if (shift >= sizeof(fixed) * 8 - 24)
  561.         return_error(gs_error_limitcheck);
  562.     if (vf < 0)
  563.         mantissa = -mantissa;
  564.     *pr = (fixed) (mantissa << shift);
  565.     } else
  566.     *pr = (shift < -24 ? fixed_0 :
  567.            vf < 0 ? -(fixed) (mantissa >> -shift) :        /* truncate */
  568.            (fixed) (mantissa >> -shift));
  569.     return 0;
  570. }
  571. int
  572. set_double2fixed_(fixed * pr, ulong /*double lo */ lo,
  573.           long /*double hi */ hi, int frac_bits)
  574. {
  575.     fixed mantissa;
  576.     int shift;
  577.  
  578.     if (!(hi & 0x7ff00000)) {
  579.     *pr = fixed_0;
  580.     return 0;
  581.     }
  582.     /* We only use 31 bits of mantissa even if sizeof(long) > 4. */
  583.     mantissa = (fixed) (((hi & 0xfffff) << 10) | (lo >> 22) | 0x40000000);
  584.     shift = ((hi >> 20) & 2047) - (1023 + 30) + frac_bits;
  585.     if (shift > 0)
  586.     return_error(gs_error_limitcheck);
  587.     *pr = (shift < -30 ? fixed_0 :
  588.        hi < 0 ? -(fixed) (mantissa >> -shift) :    /* truncate */
  589.        (fixed) (mantissa >> -shift));
  590.     return 0;
  591. }
  592. /*
  593.  * Given a fixed value x with fbits bits of fraction, set v to the mantissa
  594.  * (left-justified in 32 bits) and f to the exponent word of the
  595.  * corresponding floating-point value with mbits bits of mantissa in the
  596.  * first word.  (The exponent part of f is biased by -1, because we add the
  597.  * top 1-bit of the mantissa to it.)
  598.  */
  599. static const byte f2f_shifts[] =
  600. {4, 3, 2, 2, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0};
  601.  
  602. #define f2f_declare(v, f)\
  603.     ulong v;\
  604.     long f
  605. #define f2f(x, v, f, mbits, fbits)\
  606.     if ( x < 0 )\
  607.       f = 0xc0000000 + (29 << mbits) - ((long)fbits << mbits), v = -x;\
  608.     else\
  609.       f = 0x40000000 + (29 << mbits) - ((long)fbits << mbits), v = x;\
  610.     if ( v < 0x8000 )\
  611.       v <<= 15, f -= 15 << mbits;\
  612.     if ( v < 0x800000 )\
  613.       v <<= 8, f -= 8 << mbits;\
  614.     if ( v < 0x8000000 )\
  615.       v <<= 4, f -= 4 << mbits;\
  616.     { int shift = f2f_shifts[v >> 28];\
  617.       v <<= shift, f -= shift << mbits;\
  618.     }
  619. long
  620. fixed2float_(fixed x, int frac_bits)
  621. {
  622.     f2f_declare(v, f);
  623.  
  624.     if (x == 0)
  625.     return 0;
  626.     f2f(x, v, f, mbits_float, frac_bits);
  627.     return f + (((v >> 7) + 1) >> 1);
  628. }
  629. void
  630. set_fixed2double_(double *pd, fixed x, int frac_bits)
  631. {
  632.     f2f_declare(v, f);
  633.  
  634.     if (x == 0) {
  635.     ((long *)pd)[1 - arch_is_big_endian] = 0;
  636.     ((ulong *) pd)[arch_is_big_endian] = 0;
  637.     } else {
  638.     f2f(x, v, f, mbits_double, frac_bits);
  639.     ((long *)pd)[1 - arch_is_big_endian] = f + (v >> 11);
  640.     ((ulong *) pd)[arch_is_big_endian] = v << 21;
  641.     }
  642. }
  643.  
  644. #endif                /* USE_FPU_FIXED */
  645.  
  646. /*
  647.  * Compute A * B / C when 0 <= B < C and A * B exceeds (or might exceed)
  648.  * the capacity of a long.
  649.  * Note that this procedure takes the floor, rather than truncating
  650.  * towards zero, if A < 0.  This ensures that 0 <= R < C.
  651.  */
  652.  
  653. #define num_bits (sizeof(fixed) * 8)
  654. #define half_bits (num_bits / 2)
  655. #define half_mask ((1L << half_bits) - 1)
  656.  
  657. /*
  658.  * If doubles aren't wide enough, we lose too much precision by using double
  659.  * arithmetic: we have to use the slower, accurate fixed-point algorithm.
  660.  */
  661. #if USE_FPU_FIXED || (arch_double_mantissa_bits < arch_sizeof_long * 12)
  662.  
  663. #ifdef DEBUG
  664. struct {
  665.     long mnanb, mnab, manb, mab, mnc, mdq, mde, mds, mqh, mql;
  666. } fmq_stat;
  667. #  define mincr(x) ++fmq_stat.x
  668. #else
  669. #  define mincr(x) DO_NOTHING
  670. #endif
  671. fixed
  672. fixed_mult_quo(fixed signed_A, fixed B, fixed C)
  673. {
  674.     /* First compute A * B in double-fixed precision. */
  675.     ulong A = (signed_A < 0 ? -signed_A : signed_A);
  676.     long msw;
  677.     ulong lsw;
  678.     ulong p1;
  679.  
  680.     if (B <= half_mask) {
  681.     if (A <= half_mask) {
  682.         ulong P = A * B;
  683.         fixed Q = P / (ulong)C;
  684.  
  685.         mincr(mnanb);
  686.         /* If A < 0 and the division isn't exact, take the floor. */
  687.         return (signed_A >= 0 ? Q : Q * C == P ? -Q : ~Q /* -Q - 1 */);
  688.     }
  689.     /*
  690.      * We might still have C <= half_mask, which we can
  691.      * handle with a simpler algorithm.
  692.      */
  693.     lsw = (A & half_mask) * B;
  694.     p1 = (A >> half_bits) * B;
  695.     if (C <= half_mask) {
  696.         fixed q0 = (p1 += lsw >> half_bits) / C;
  697.         ulong rem = ((p1 - C * q0) << half_bits) + (lsw & half_mask);
  698.         ulong q1 = rem / (ulong)C;
  699.         fixed Q = (q0 << half_bits) + q1;
  700.  
  701.         mincr(mnc);
  702.         /* If A < 0 and the division isn't exact, take the floor. */
  703.         return (signed_A >= 0 ? Q : q1 * C == rem ? -Q : ~Q);
  704.     }
  705.     msw = p1 >> half_bits;
  706.     mincr(manb);
  707.     } else if (A <= half_mask) {
  708.     p1 = A * (B >> half_bits);
  709.     msw = p1 >> half_bits;
  710.     lsw = A * (B & half_mask);
  711.     mincr(mnab);
  712.     } else {            /* We have to compute all 4 products.  :-( */
  713.     ulong lo_A = A & half_mask;
  714.     ulong hi_A = A >> half_bits;
  715.     ulong lo_B = B & half_mask;
  716.     ulong hi_B = B >> half_bits;
  717.     ulong p1x = hi_A * lo_B;
  718.  
  719.     msw = hi_A * hi_B;
  720.     lsw = lo_A * lo_B;
  721.     p1 = lo_A * hi_B;
  722.     if (p1 > max_ulong - p1x)
  723.         msw += 1L << half_bits;
  724.     p1 += p1x;
  725.     msw += p1 >> half_bits;
  726.     mincr(mab);
  727.     }
  728.     /* Finish up by adding the low half of p1 to the high half of lsw. */
  729. #if max_fixed < max_long
  730.     p1 &= half_mask;
  731. #endif
  732.     p1 <<= half_bits;
  733.     if (p1 > max_ulong - lsw)
  734.     msw++;
  735.     lsw += p1;
  736.     /*
  737.      * Now divide the double-length product by C.  Note that we know msw
  738.      * < C (otherwise the quotient would overflow).  Start by shifting
  739.      * (msw,lsw) and C left until C >= 1 << (num_bits - 1).
  740.      */
  741.     {
  742.     ulong denom = C;
  743.     int shift = 0;
  744.  
  745. #define bits_4th (num_bits / 4)
  746.     if (denom < 1L << (num_bits - bits_4th)) {
  747.         mincr(mdq);
  748.         denom <<= bits_4th, shift += bits_4th;
  749.     }
  750. #undef bits_4th
  751. #define bits_8th (num_bits / 8)
  752.     if (denom < 1L << (num_bits - bits_8th)) {
  753.         mincr(mde);
  754.         denom <<= bits_8th, shift += bits_8th;
  755.     }
  756. #undef bits_8th
  757.     while (!(denom & (-1L << (num_bits - 1)))) {
  758.         mincr(mds);
  759.         denom <<= 1, ++shift;
  760.     }
  761.     msw = (msw << shift) + (lsw >> (num_bits - shift));
  762.     lsw <<= shift;
  763. #if max_fixed < max_long
  764.     lsw &= (1L << (sizeof(fixed) * 8)) - 1;
  765. #endif
  766.     /* Compute a trial upper-half quotient. */
  767.     {
  768.         ulong hi_D = denom >> half_bits;
  769.         ulong lo_D = denom & half_mask;
  770.         ulong hi_Q = (ulong) msw / hi_D;
  771.  
  772.         /* hi_Q might be too high by 1 or 2, but it isn't too low. */
  773.         ulong p0 = hi_Q * hi_D;
  774.         ulong p1 = hi_Q * lo_D;
  775.         ulong hi_P;
  776.  
  777.         while ((hi_P = p0 + (p1 >> half_bits)) > msw ||
  778.            (hi_P == msw && ((p1 & half_mask) << half_bits) > lsw)
  779.         ) {        /* hi_Q was too high by 1. */
  780.         --hi_Q;
  781.         p0 -= hi_D;
  782.         p1 -= lo_D;
  783.         mincr(mqh);
  784.         }
  785.         p1 = (p1 & half_mask) << half_bits;
  786.         if (p1 > lsw)
  787.         msw--;
  788.         lsw -= p1;
  789.         msw -= hi_P;
  790.         /* Now repeat to get the lower-half quotient. */
  791.         msw = (msw << half_bits) + (lsw >> half_bits);
  792. #if max_fixed < max_long
  793.         lsw &= half_mask;
  794. #endif
  795.         lsw <<= half_bits;
  796.         {
  797.         ulong lo_Q = (ulong) msw / hi_D;
  798.         long Q;
  799.  
  800.         p1 = lo_Q * lo_D;
  801.         p0 = lo_Q * hi_D;
  802.         while ((hi_P = p0 + (p1 >> half_bits)) > msw ||
  803.                (hi_P == msw && ((p1 & half_mask) << half_bits) > lsw)
  804.             ) {        /* lo_Q was too high by 1. */
  805.             --lo_Q;
  806.             p0 -= hi_D;
  807.             p1 -= lo_D;
  808.             mincr(mql);
  809.         }
  810.         Q = (hi_Q << half_bits) + lo_Q;
  811.         return (signed_A >= 0 ? Q : p0 | p1 ? ~Q /* -Q - 1 */ : -Q);
  812.         }
  813.     }
  814.     }
  815. }
  816.  
  817. #else                /* can approximate using doubles */
  818.  
  819. /*
  820.  * Compute A * B / C as above.  Since a double doesn't have enough bits to
  821.  * represent the product of two longs, we have to do it in two steps.
  822.  */
  823. fixed
  824. fixed_mult_quo(fixed signed_A, fixed B, fixed C)
  825. {
  826. #define MAX_OTHER_FACTOR\
  827.   (1L << (arch_double_mantissa_bits - sizeof(fixed) * 8))
  828.     if (B < MAX_OTHER_FACTOR || any_abs(signed_A) < MAX_OTHER_FACTOR) {
  829.     /* The double computation will be exact. */
  830.     return (fixed)floor((double)signed_A * B / C);
  831.     }
  832. #undef MAX_OTHER_FACTOR
  833.     {
  834.     /* Use 2 double steps. */
  835.     fixed bhi = B >> half_bits;
  836.     fixed qhi = (fixed)floor((double)signed_A * bhi / C);
  837.     fixed rhi = signed_A * bhi - qhi * C;
  838.     fixed blo = B & half_mask;
  839.     fixed qlo =
  840.         (fixed)floor(((double)rhi * (1L << half_bits) +
  841.               (double)signed_A * blo) / C);
  842.  
  843.     return (qhi << half_bits) + qlo;
  844.     }
  845. }
  846.  
  847. #endif
  848.  
  849. #undef num_bits
  850. #undef half_bits
  851. #undef half_mask
  852.  
  853. /* Trace calls on sqrt when debugging. */
  854. double
  855. gs_sqrt(double x, const char *file, int line)
  856. {
  857.     if (gs_debug_c('~')) {
  858.     dprintf3("[~]sqrt(%g) at %s:%d\n", x, (const char *)file, line);
  859.     dflush();
  860.     }
  861.     return orig_sqrt(x);
  862. }
  863.  
  864. /*
  865.  * Define sine and cosine functions that take angles in degrees rather than
  866.  * radians, and that are implemented efficiently on machines with slow
  867.  * (or no) floating point.
  868.  */
  869. #if USE_FPU < 0            /****** maybe should be <= 0 ? ***** */
  870.  
  871. #define sin0 0.00000000000000000
  872. #define sin1 0.01745240643728351
  873. #define sin2 0.03489949670250097
  874. #define sin3 0.05233595624294383
  875. #define sin4 0.06975647374412530
  876. #define sin5 0.08715574274765817
  877. #define sin6 0.10452846326765346
  878. #define sin7 0.12186934340514748
  879. #define sin8 0.13917310096006544
  880. #define sin9 0.15643446504023087
  881. #define sin10 0.17364817766693033
  882. #define sin11 0.19080899537654480
  883. #define sin12 0.20791169081775931
  884. #define sin13 0.22495105434386498
  885. #define sin14 0.24192189559966773
  886. #define sin15 0.25881904510252074
  887. #define sin16 0.27563735581699916
  888. #define sin17 0.29237170472273671
  889. #define sin18 0.30901699437494740
  890. #define sin19 0.32556815445715670
  891. #define sin20 0.34202014332566871
  892. #define sin21 0.35836794954530027
  893. #define sin22 0.37460659341591201
  894. #define sin23 0.39073112848927377
  895. #define sin24 0.40673664307580015
  896. #define sin25 0.42261826174069944
  897. #define sin26 0.43837114678907740
  898. #define sin27 0.45399049973954675
  899. #define sin28 0.46947156278589081
  900. #define sin29 0.48480962024633706
  901. #define sin30 0.50000000000000000
  902. #define sin31 0.51503807491005416
  903. #define sin32 0.52991926423320490
  904. #define sin33 0.54463903501502708
  905. #define sin34 0.55919290347074679
  906. #define sin35 0.57357643635104605
  907. #define sin36 0.58778525229247314
  908. #define sin37 0.60181502315204827
  909. #define sin38 0.61566147532565829
  910. #define sin39 0.62932039104983739
  911. #define sin40 0.64278760968653925
  912. #define sin41 0.65605902899050728
  913. #define sin42 0.66913060635885824
  914. #define sin43 0.68199836006249848
  915. #define sin44 0.69465837045899725
  916. #define sin45 0.70710678118654746
  917. #define sin46 0.71933980033865108
  918. #define sin47 0.73135370161917046
  919. #define sin48 0.74314482547739413
  920. #define sin49 0.75470958022277201
  921. #define sin50 0.76604444311897801
  922. #define sin51 0.77714596145697090
  923. #define sin52 0.78801075360672190
  924. #define sin53 0.79863551004729283
  925. #define sin54 0.80901699437494745
  926. #define sin55 0.81915204428899180
  927. #define sin56 0.82903757255504174
  928. #define sin57 0.83867056794542394
  929. #define sin58 0.84804809615642596
  930. #define sin59 0.85716730070211222
  931. #define sin60 0.86602540378443860
  932. #define sin61 0.87461970713939574
  933. #define sin62 0.88294759285892688
  934. #define sin63 0.89100652418836779
  935. #define sin64 0.89879404629916704
  936. #define sin65 0.90630778703664994
  937. #define sin66 0.91354545764260087
  938. #define sin67 0.92050485345244037
  939. #define sin68 0.92718385456678731
  940. #define sin69 0.93358042649720174
  941. #define sin70 0.93969262078590832
  942. #define sin71 0.94551857559931674
  943. #define sin72 0.95105651629515353
  944. #define sin73 0.95630475596303544
  945. #define sin74 0.96126169593831889
  946. #define sin75 0.96592582628906831
  947. #define sin76 0.97029572627599647
  948. #define sin77 0.97437006478523525
  949. #define sin78 0.97814760073380558
  950. #define sin79 0.98162718344766398
  951. #define sin80 0.98480775301220802
  952. #define sin81 0.98768834059513777
  953. #define sin82 0.99026806874157036
  954. #define sin83 0.99254615164132198
  955. #define sin84 0.99452189536827329
  956. #define sin85 0.99619469809174555
  957. #define sin86 0.99756405025982420
  958. #define sin87 0.99862953475457383
  959. #define sin88 0.99939082701909576
  960. #define sin89 0.99984769515639127
  961. #define sin90 1.00000000000000000
  962.  
  963. private const double sin_table[361] =
  964. {
  965.     sin0,
  966.     sin1, sin2, sin3, sin4, sin5, sin6, sin7, sin8, sin9, sin10,
  967.     sin11, sin12, sin13, sin14, sin15, sin16, sin17, sin18, sin19, sin20,
  968.     sin21, sin22, sin23, sin24, sin25, sin26, sin27, sin28, sin29, sin30,
  969.     sin31, sin32, sin33, sin34, sin35, sin36, sin37, sin38, sin39, sin40,
  970.     sin41, sin42, sin43, sin44, sin45, sin46, sin47, sin48, sin49, sin50,
  971.     sin51, sin52, sin53, sin54, sin55, sin56, sin57, sin58, sin59, sin60,
  972.     sin61, sin62, sin63, sin64, sin65, sin66, sin67, sin68, sin69, sin70,
  973.     sin71, sin72, sin73, sin74, sin75, sin76, sin77, sin78, sin79, sin80,
  974.     sin81, sin82, sin83, sin84, sin85, sin86, sin87, sin88, sin89, sin90,
  975.     sin89, sin88, sin87, sin86, sin85, sin84, sin83, sin82, sin81, sin80,
  976.     sin79, sin78, sin77, sin76, sin75, sin74, sin73, sin72, sin71, sin70,
  977.     sin69, sin68, sin67, sin66, sin65, sin64, sin63, sin62, sin61, sin60,
  978.     sin59, sin58, sin57, sin56, sin55, sin54, sin53, sin52, sin51, sin50,
  979.     sin49, sin48, sin47, sin46, sin45, sin44, sin43, sin42, sin41, sin40,
  980.     sin39, sin38, sin37, sin36, sin35, sin34, sin33, sin32, sin31, sin30,
  981.     sin29, sin28, sin27, sin26, sin25, sin24, sin23, sin22, sin21, sin20,
  982.     sin19, sin18, sin17, sin16, sin15, sin14, sin13, sin12, sin11, sin10,
  983.     sin9, sin8, sin7, sin6, sin5, sin4, sin3, sin2, sin1, sin0,
  984.     -sin1, -sin2, -sin3, -sin4, -sin5, -sin6, -sin7, -sin8, -sin9, -sin10,
  985.     -sin11, -sin12, -sin13, -sin14, -sin15, -sin16, -sin17, -sin18, -sin19, -sin20,
  986.     -sin21, -sin22, -sin23, -sin24, -sin25, -sin26, -sin27, -sin28, -sin29, -sin30,
  987.     -sin31, -sin32, -sin33, -sin34, -sin35, -sin36, -sin37, -sin38, -sin39, -sin40,
  988.     -sin41, -sin42, -sin43, -sin44, -sin45, -sin46, -sin47, -sin48, -sin49, -sin50,
  989.     -sin51, -sin52, -sin53, -sin54, -sin55, -sin56, -sin57, -sin58, -sin59, -sin60,
  990.     -sin61, -sin62, -sin63, -sin64, -sin65, -sin66, -sin67, -sin68, -sin69, -sin70,
  991.     -sin71, -sin72, -sin73, -sin74, -sin75, -sin76, -sin77, -sin78, -sin79, -sin80,
  992.     -sin81, -sin82, -sin83, -sin84, -sin85, -sin86, -sin87, -sin88, -sin89, -sin90,
  993.     -sin89, -sin88, -sin87, -sin86, -sin85, -sin84, -sin83, -sin82, -sin81, -sin80,
  994.     -sin79, -sin78, -sin77, -sin76, -sin75, -sin74, -sin73, -sin72, -sin71, -sin70,
  995.     -sin69, -sin68, -sin67, -sin66, -sin65, -sin64, -sin63, -sin62, -sin61, -sin60,
  996.     -sin59, -sin58, -sin57, -sin56, -sin55, -sin54, -sin53, -sin52, -sin51, -sin50,
  997.     -sin49, -sin48, -sin47, -sin46, -sin45, -sin44, -sin43, -sin42, -sin41, -sin40,
  998.     -sin39, -sin38, -sin37, -sin36, -sin35, -sin34, -sin33, -sin32, -sin31, -sin30,
  999.     -sin29, -sin28, -sin27, -sin26, -sin25, -sin24, -sin23, -sin22, -sin21, -sin20,
  1000.     -sin19, -sin18, -sin17, -sin16, -sin15, -sin14, -sin13, -sin12, -sin11, -sin10,
  1001.     -sin9, -sin8, -sin7, -sin6, -sin5, -sin4, -sin3, -sin2, -sin1, -sin0
  1002. };
  1003.  
  1004. double
  1005. gs_sin_degrees(double ang)
  1006. {
  1007.     int ipart;
  1008.  
  1009.     if (is_fneg(ang))
  1010.     ang = 180 - ang;
  1011.     ipart = (int)ang;
  1012.     if (ipart >= 360) {
  1013.     int arem = ipart % 360;
  1014.  
  1015.     ang -= (ipart - arem);
  1016.     ipart = arem;
  1017.     }
  1018.     return
  1019.     (ang == ipart ? sin_table[ipart] :
  1020.      sin_table[ipart] + (sin_table[ipart + 1] - sin_table[ipart]) *
  1021.      (ang - ipart));
  1022. }
  1023.  
  1024. double
  1025. gs_cos_degrees(double ang)
  1026. {
  1027.     int ipart;
  1028.  
  1029.     if (is_fneg(ang))
  1030.     ang = 90 - ang;
  1031.     else
  1032.     ang += 90;
  1033.     ipart = (int)ang;
  1034.     if (ipart >= 360) {
  1035.     int arem = ipart % 360;
  1036.  
  1037.     ang -= (ipart - arem);
  1038.     ipart = arem;
  1039.     }
  1040.     return
  1041.     (ang == ipart ? sin_table[ipart] :
  1042.      sin_table[ipart] + (sin_table[ipart + 1] - sin_table[ipart]) *
  1043.      (ang - ipart));
  1044. }
  1045.  
  1046. void
  1047. gs_sincos_degrees(double ang, gs_sincos_t * psincos)
  1048. {
  1049.     psincos->sin = gs_sin_degrees(ang);
  1050.     psincos->cos = gs_cos_degrees(ang);
  1051.     psincos->orthogonal =
  1052.     (is_fzero(psincos->sin) || is_fzero(psincos->cos));
  1053. }
  1054.  
  1055. #else /* we have floating point */
  1056.  
  1057. static const int isincos[5] =
  1058. {0, 1, 0, -1, 0};
  1059.  
  1060. double
  1061. gs_sin_degrees(double ang)
  1062. {
  1063.     double quot = ang / 90;
  1064.  
  1065.     if (floor(quot) == quot) {
  1066.     /*
  1067.      * We need 4.0, rather than 4, here because of non-ANSI compilers.
  1068.      * The & 3 is because quot might be negative.
  1069.      */
  1070.     return isincos[(int)fmod(quot, 4.0) & 3];
  1071.     }
  1072.     return sin(ang * (M_PI / 180));
  1073. }
  1074.  
  1075. double
  1076. gs_cos_degrees(double ang)
  1077. {
  1078.     double quot = ang / 90;
  1079.  
  1080.     if (floor(quot) == quot) {
  1081.     /* See above re the following line. */
  1082.     return isincos[((int)fmod(quot, 4.0) & 3) + 1];
  1083.     }
  1084.     return cos(ang * (M_PI / 180));
  1085. }
  1086.  
  1087. void
  1088. gs_sincos_degrees(double ang, gs_sincos_t * psincos)
  1089. {
  1090.     double quot = ang / 90;
  1091.  
  1092.     if (floor(quot) == quot) {
  1093.     /* See above re the following line. */
  1094.     int quads = (int)fmod(quot, 4.0) & 3;
  1095.  
  1096.     psincos->sin = isincos[quads];
  1097.     psincos->cos = isincos[quads + 1];
  1098.     psincos->orthogonal = true;
  1099.     } else {
  1100.     double arad = ang * (M_PI / 180);
  1101.  
  1102.     psincos->sin = sin(arad);
  1103.     psincos->cos = cos(arad);
  1104.     psincos->orthogonal = false;
  1105.     }
  1106. }
  1107.  
  1108. #endif /* USE_FPU */
  1109.  
  1110. /*
  1111.  * Define an atan2 function that returns an angle in degrees and uses
  1112.  * the PostScript quadrant rules.  Note that it may return
  1113.  * gs_error_undefinedresult.
  1114.  */
  1115. int
  1116. gs_atan2_degrees(double y, double x, double *pangle)
  1117. {
  1118.     if (y == 0) {    /* on X-axis, special case */
  1119.     if (x == 0)
  1120.         return_error(gs_error_undefinedresult);
  1121.     *pangle = (x < 0 ? 180 : 0);
  1122.     } else {
  1123.     double result = atan2(y, x) * radians_to_degrees;
  1124.  
  1125.     if (result < 0)
  1126.         result += 360;
  1127.     *pangle = result;
  1128.     }
  1129.     return 0;
  1130. }
  1131.